%%%%%%%%%%%%%%%%%%
% option return factor model+ some stats
% model: Return = Structural factors (HC, VR, JR)+[HRP, VRP] 
% Liuren Wu, liuren.wu@baruch.cuny.edu
%%%%%%%%%%%%%%%%%%
function [Bi,nobt,Af,Bf,R2s,stats,SXV,SXV2,bi,R2i,qx,qy,qv,VB] =FunDailyORSestimation(t,y,x,HVV,IV,ids,dds,dd,nlag,nx,winth,minN)

tt=dds==dd(t);  sdt=ids(tt);
yt=y(tt);

Nt=sum(tt);
rsv=NaN(Nt,nlag);
jj=1:min(t-1,nlag);
nj=length(jj);
Biv=NaN(3,nj);Bi0=ones(3,1)/3;
for j=1:nj
    ll=jj(j);
    ttp=dds==dd(t-ll);sdp=ids(ttp);ytp=y(ttp);
    [~,ind1,ind2]=intersect(sdt,sdp);
    rsv(ind1,j)=ytp(ind2);
    
    HVVp=HVV(ttp,:);
    
    tti=dds==dd(t-ll+1);HVi=HVV(tti,1);sdi=ids(tti);
    [~,ind1,ind2]=intersect(sdi,sdp);
    yi=HVi(ind1);
    Xi=HVVp(ind2,:);
    ii=isfinite(yi+sum(Xi,2)); yi=yi(ii);Xi=Xi(ii,:);
    XX=Xi'*Xi;
    if rank(XX)<3
        Pi=diag(diag(XX))*1e-4;
        Biv(:,j)=(XX+Pi)\(Xi'*yi + Pi*Bi0);
    else
        Biv(:,j)=Xi\yi;
    end
end
Bi=nanmean(Biv,2);

HVVt=HVV(tt,:);FHV=HVVt*Bi;
VRP=100*(IV(tt) -FHV); % VRP
HRP=nanmean(rsv(:,2:nlag),2);%HRP

tvrx=(nanstd(rsv,0,2)); %rolling-window return volatility estimator

xt=[x(tt,:),HRP,VRP]; %xt(:,3)=xt(:,3).*FHV;

nobt=[Nt,sum(isfinite(yt)),sum(isfinite(xt))];
xt=winsorize(xt,winth);
tvr=winsorize(tvrx,winth);
%1. Descriptor stats
Zt=[yt,xt,tvr];
stats=[...
    nanmean(Zt);
    nanstd(Zt,0);
    prctile(Zt,[10,25,50,75,90]');
    skewness(Zt,0,1);
    kurtosis(Zt,0,1)-3;
    100*corr(yt,Zt,'rows','pairwise');
    100*corr(yt,Zt,'rows','pairwise','type','Spearman')];
%quintiles
dv=[yt,xt]*NaN;nmd=5;
for j=1:nx+1
    if sum(isfinite(Zt(:,j)))>minN
    dv(:,j) = discretize(Zt(:,j),quantile(Zt(:,j),[0:nmd]/nmd));
    end
end

qx=NaN(nx,nmd);qy=qx;qv=qy;
for k=1:nmd
    qx(:,k)=nanmean(xt(dv(:,1)==k,:));
    for j=1:nx
        qy(j,k)=nanmean(yt(dv(:,j+1)==k));
        qv(j,k)=nanmean(tvrx(dv(:,j+1)==k));
    end
end

ii=isfinite(sum(xt(:,1:3),2)+tvr);
yi=standardize(tvr(ii),1);
xi=standardize(xt(ii,1:3),1);
bi=xi\yi;
ei=yi-xi*bi;
R2i=1-var(ei);


% %regression
SZt=NaN(Nt,nx); 
indfx=zeros(nx,1);
SXV=NaN(1,nx); 
SZ2=SZt;
for j=1:nx
    indf=isfinite(xt(:,j)+yt); Nj=sum(indf);
    if Nj>minN
        indfx(j)=1;
        xtj=xt(indf,j);sxj=std(xtj);mxj=mean(xtj);
        if j<=3
            SZt(indf,j)=xtj./sxj;
        else
            SZt(indf,j)=(xtj-mxj)./sxj;
        end
        SXV(j)=sxj;
        SZ2(indf,j)=xtj;
    end
end


Bf=NaN(nx,1); R2s=NaN(1,1);Af=NaN(1,1); VB=Bf;
indff=find(indfx==1);
if ~isempty(indff)
    Ft = SZt(:,indff);     [Af(1),Bf(indff,1),R2s,VB(indff),SXV2]=FunFMRegression(Ft,yt,SZ2);
end

end

function [A,B,AR2,vb,CO]=FunFMRegression(Ft,yt,SZ2)
indf=isfinite(sum(Ft,2)+yt); Nf=sum(indf);

Xtf=[ones(Nf,1),Ft(indf,:)];ytf=yt(indf); 
B=Xtf\ytf; 
yh=Xtf*B; 
ef=ytf-yh;
nx=size(Xtf,2);
AR2= 1-(sum(ef.^2)./(Nf-nx))/var(ytf);

A=B(1); B=B(2:end);
Stf=SZ2(indf,:);
C=Xtf\Stf;
CO=diag(C(2:end,:));

vb=B.*(cov(Xtf(:,2:end))*B)./var(ytf);



end




